%% thu nghiem 2:
% tinh xac suat x co trong gallery roi tinh AP va MAP
% trong phan nay ta su dung du lieu anh fa, ql va ket qua cua viec training
% fa va ql o trang thai tot nhat ung voi factor 32
function MAP = FaceRetrieval_ThuNghiem02
    % sinh nghia cac mang luu gia tri tuong ung.
    dieuKien = {'fa','ql'};
    
    % bien giua gia tri cua tat ca cac x. Nhu vay sau dong for ben duoi thi
    % ket qua c?a allpixelData se la 14700*638 (tuc la 319 x 2) do ta xem
    % fa va ql coi nhu la 2 tap khac biet chua 638 nguoi khac nhau.
    allpixelData=[];
    
    %% duyet tung dieu kien load thong tin tuong ung vao mang allpixelData
    for(dk = 1:length(dieuKien))
        fileName = sprintf('./PixelDataTri_%s.mat',dieuKien{dk});
        load(fileName);
        allpixelData=[allpixelData pixelData];
    end
    nImage = size(allpixelData,2);
    %% Tao danh sach cac X.
    for(Xi =1:size(allpixelData,2)) % 638 ca nhan
        % moi mot ca nhan co 1 vector dac trung goi la X
        X{Xi} = allpixelData(:,Xi);
    end
    
    %% tinh xac suat  1 xi so voi tat ca cac xj con lai
    % Ta se dung lai mo hinh xac suat cua tac gia: Simon J.D. Prince. Trong
    % paper: Tied Factor Analysis for Face Recognition across Large Pose Differences
    % 
    fprintf('Dang su dung mo hinh 32 factor!\n');
    
    % load file tuong ung de lay F, sigma va meanVector
    outName = 'faql';
    filename = sprintf('./TiedFacModelPix_%s_%d.mat',outName,32);
    load(filename,'meanVec','FEst','sigmaEst','FEstTied','sigmaEstTied','meanVecTied');
    
    %% tinh lai sigma giong nhu tac gia
    sigmaEstTied{1} = sigmaEstTied{1}*32*10/48;
    sigmaEstTied{2} = sigmaEstTied{2}*32*10/48;
    
    %% su dung lai mo hinh tinh xac suat cua tac gia.
    for(irun=1:nImage)
        % tinh lai tap gellary = allpixelData - X{irun}
        fprintf('dang tinh xac suat cho anh thu: %d\n',irun);
        if(irun == 1)
            gellaryData = allpixelData(:,irun+1:nImage);
        elseif (irun == nImage)
            gellaryData = allpixelData(:,1:irun-1);
        else
            gellaryData = allpixelData(:,1:irun-1);
            gellaryData =[gellaryData allpixelData(:,irun+1:nImage)];
        end
        modelLogLikelihoods = getLogLikeIDTied(FEstTied,sigmaEstTied,meanVecTied,gellaryData,X{irun});
        postProb = modelLogLikelihoods-repmat(max(modelLogLikelihoods),nImage-1,1); % nImage-1: vi phai loai tu the dang xet
        postProb = exp(postProb);
        postProb = postProb./repmat(sum(postProb),nImage-1,1);
        Pro{irun} = postProb;
    end
    
    % sau khi tinh duoc xac suat thi ta tien hanh tim vi tri AP
    for(irun=1:nImage)
        giaTri{irun} = Pro{irun}((nImage/2)+(irun-1));
        
    end

    % sap xep cho xac xuat giam dan
    for(irun=1:nImage)
        Pro{irun} = sort(Pro{irun},'descend');
    end
    
    % tim gia tri AP
    for(irun=1:nImage)
        vitri = find(Pro{1}== giaTri{irun});
        AP{irun} = 1/vitri;
    end
    
    %% tinh trung binh ap
    MAP = 0;
    for(irun = 1:lenght(AP))
        MAP = MAP + AP{irun};
    end
    MAP = MAP/length(AP);
    %% ***************************************************************************************
function modelLogLikelihoods = getLogLikeIDTied(FEst,sigmaEst,meanVec,galleryData,probeData);

 
%define space for log likelihoods
[nDim nImage] = size(galleryData);
modelLogLikelihoods = zeros(nImage,1);

%remove means from data
galleryData = galleryData-repmat(meanVec{1},[1,nImage]);
probeData = probeData-repmat(meanVec{2},[1,1]);%%

%preprocess factor analysis model
HIGHEST_N = 2; %highest number of faces explained by one variable
factorModelTied = preProcessFATiedPair(FEst,sigmaEst);
galleryDataPP = preProcessFactorDataTied(factorModelTied,galleryData,1);
probeDataPP = preProcessFactorDataTied(factorModelTied,probeData,2);

%precalculate single probability terms
logLikeGalleryNoMatch = zeros(nImage,1);
for (cImage = 1:nImage)
    logLikeGalleryNoMatch(cImage) = getLogLikeMatchTiedFA(factorModelTied,galleryDataPP(cImage),1);
end;

%for each probe
%for (cProbe = 1:nImage);
    %for each gallery image
    cProbe = 1;
    for(cGallery = 1:nImage);
        %calculate probability of this model
        %fprintf('cp=%d,cg=%d\n',cProbe,cGallery);
        logLikeModel= sum(logLikeGalleryNoMatch(:))-logLikeGalleryNoMatch(cGallery);
        logLikeModel = logLikeModel+getLogLikeMatchTiedFA(factorModelTied,[galleryDataPP(cGallery) probeDataPP(cProbe)],1:2);
        modelLogLikelihoods(cGallery,cProbe) = logLikeModel;
    end;
%end;
drawnow;



%==========================================================================
%==========================================================================

function factorModel = preProcessFATiedPair(FAll,covDiagAll);

[N_DATA_DIM N_HIDDEN_DIM] = size(FAll{1});

factorModel.F = FAll;
nCond = size(FAll,2);
for (cCond = 1:nCond)
    factorModel.F{cCond} = FAll{cCond};
    factorModel.covDiag{cCond} = covDiagAll{cCond};
    factorModel.invCovDiag{cCond} = 1./covDiagAll{cCond};
    factorModel.FWeighted{cCond} = (FAll{cCond}.*repmat(factorModel.invCovDiag{cCond},1,N_HIDDEN_DIM));
    factorModel.invFSFPlusIDiag{cCond} = inv(factorModel.F{cCond}'*factorModel.FWeighted{cCond}+eye(N_HIDDEN_DIM));
    factorModel.detInvFSFPlusIDiag{cCond} = det(factorModel.invFSFPlusIDiag{cCond});
end;

factorModel.invSumFSFPlusIDiag = inv(factorModel.F{1}'*factorModel.FWeighted{1} ...
                +factorModel.F{2}'*factorModel.FWeighted{2}+eye(N_HIDDEN_DIM));

%factorModel.invSumFSFPlusIDiag = inv(2*factorModel.F{1}'*factorModel.FWeighted{1}+eye(N_HIDDEN_DIM));
            
factorModel.detInvSumFSFPlusIDiag = det(factorModel.invSumFSFPlusIDiag);


%==========================================================================
%==========================================================================
function dataPP = preProcessFactorDataTied(factorModel,data,condition);

%condition = 1;
[N_DATA_DIM N_DATA] = size(data);
for (cData = 1:N_DATA);
    dataPP(cData).FTinvSx = factorModel.F{condition}'*(data(:,cData).*factorModel.invCovDiag{condition});
    dataPP(cData).logProbTerm = getLogGaussProbability(data(:,cData),0,factorModel.covDiag{condition});
end;
    
%==========================================================================
%==========================================================================
function logLike = getLogLikeMatchTiedFA(factorModel,data, conditions);

nData = length(data);
N_HIDDEN = size(factorModel.F{1},2);
%conditions = ones(size(conditions));
if (length(conditions)==1)
    logLike = 0.5*log(factorModel.detInvFSFPlusIDiag{conditions});
    sumWeightedData = zeros(N_HIDDEN,1);

    logLike = logLike+data.logProbTerm;
    sumWeighted = data.FTinvSx;
    logLike = logLike+0.5*(sumWeighted'*factorModel.invFSFPlusIDiag{conditions}*sumWeighted);
    
elseif (length(conditions==2))
    logLike = 0.5*log(factorModel.detInvSumFSFPlusIDiag);
    logLike = logLike+data(1).logProbTerm;
    logLike = logLike+data(2).logProbTerm;
    sumWeighted = data(1).FTinvSx+data(2).FTinvSx;
    logLike = logLike+0.5*(sumWeighted'*factorModel.invSumFSFPlusIDiag*sumWeighted);

    drawnow;
else
    fprintf('Not written general likelihood comp yet..');
end;
 %==========================================================================
%==========================================================================

function soln = getLogGaussProbability(data,mean,covar)
%return the gaussian probability 

% The dimension of the data
dataDim = size(data,1);
%subtract the mean from the data
data = data-mean;


%if covar is diagonal then easier
if (size(covar,2)==1)
    % The result
    soln = -(dataDim/2)*log(2*pi) - (1/2)*sum(log(covar)) - ...
    0.5*data'*(data.*(1./covar));
else
    %full covariance matrix
    soln = -(dataDim/2)*log(2*pi) - (1/2)*log(det(covar)) - ...
    0.5*data'*inv(covar)*data;
end;
